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Abstract 

We study cubature formulas for d-dimensional integrals with an arbitrary 
symmetric weight function of product form. We present a construction that 
yields a high polynomial exactness: for fixed degree t = 5 or t = 7 and large 
dimension d the number of knots is only slightly larger than the lower bound 
of Moller and much smaller compared to the known constructions. 

We also show, for any odd degree I = 2k + 1, that the minimal number of 
points is almost independent of the weight function. This is also true for the 
integration over the (Euclidean) sphere. 
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1 Introduction 

Let us start with a special case of our results: We find cubature formulas with 

N(5,d,l) = d 2 + 7d+l, and N{7, d, 1) = (d 3 + 2ld 2 + 20d + 3)/3 
points such that the integral 

W)= I /(x)rfx 
J[-i,i] d 

is exactly computed for all polynomials of degree at most 5 or 7, respectively. This 
improves the known cubature formulas for degree 5 and d > 8 and for degree 7 with 
d > 10. The lower bound of Moller (1979) takes the form 

(1) iV min (5, d, 1) > d 2 + d + 1 and iV min (7, d, 1) > (d 3 + 3d 2 + 8d)/3. 
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Hence, for our method, we obtain 

(2) iV(5, d, 1) « iV min (5, d, 1) and iV(7, d, 1) « iV min (7, d, 1). 
We use ~ to denote the strong equivalence of sequences, i.e., 

v n « w n iff lim t/ n /ro n = 1. 

n— >oo 

The best results (for large <i) from the literature, see Stroud (1971) and the online 
tables of Cools, see Cools (2003), are given by 

(3) N old (5,d, 1) = 2d 2 + l and iV old (7, d, 1) = (Ad 3 - 6d 2 + Ud + 3)/3. 
More generally, we study cubature formulas 

n 

(4) Q n (f) = J2ai /(xi), o< e R, Xj G O, 

i=l 

for d-dimensional integrals 

(5) /«(/)= / /(x) ff (x)dx. 

Concerning the integral we always assume 

fl = Oi x ■ ■ ■ x fl,i 

with symmetric (and possibly unbounded) intervals Qj C 1R and the product form 

£>(x) = £i(xi) . . . Q d {x d ) 

of the weight function g. We assume that the gt are symmetric, 

g { (x) = gi(-x) 

with gi > and integrability of all polynomials, although these assumptions can be 
relaxed. Some of our results can be slightly improved in the fully symmetric case 
where, in addition, all the g^ coincide. 

Let F(£, d) be the space of all polynomials in d variables of (total) degree at 
most I. A cubature formula Q n has a degree £ of exactness if 

Qn(f) = l g d (f), V/eP(M). 

We define 

N min (£,d, g) 

to be the minimal number n of knots needed by any cubature formula Q n of degree 
£ of exactness. 

The numbers N m i n (£, d, g) and corresponding cubature formulas are only known 
in exceptional cases, see, e.g., Schmid (1983), Berens, Schmid, Xu (1995), and Cools 
(1997). Thus one is interested in upper and lower bounds for this quantity. 

One is often interested in cubature formulas with knots inside the domain and 
positive weights. While Xj G O can always be satisfied by our method, we usually 
have positive and negative weights. Actually we request Xj 6 Q, see (jH), although 
the lower of Moller also holds without this assumption. 
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2 Problem, Main Results, and Conjecture 



The lower bound of Moller (1979) for centrally symmetric weight functions is the 
following: If k is odd then 



N min (2k + l,d, ^) > 2 dim P e (M) = ^ + ^ (* + * 

If k is even then 

iV min (2A; + l,d,g)>2 dimP D (fc, d) - 1 = + ^ + ^(1 - T~ d ) ( 



s + k — 1 

s 



Here F e (k,d) denotes the subspace of P(/c, d) generated by even polynomials and 
P (/c, d) is the subspace generated by odd polynomials. We obtain (JTJ) as special 
cases and for large d the lower bounds are of the order 

2d k 

k\ ' 

See the book Mysovskikh (1981) or Cools (1997) and, for the explicit formula, Lu, 
Darmofal (2004). 

The best upper bounds were of the form 

(6) „ **. 

They can be proved with "fully symmetric formulas" (if the Qi are equal) or (in the 
general case) with the "Smolyak method" or with "sparse grids" . All these notions 
are very much related, see Section 01 Even for special weight functions g and/or for 
special I = 2k + 1 better bounds were not known. Hence there is a gap between the 
lower and the upper bound of a factor of 2 k ~ 1 and we only knew (before we wrote 
this paper) of one exception: For the weight function 

(7) £(x) = exp(— ||x|||) 



it is known for i = 5 that 

(8) d 2 + 3d + 3 

function values are enough, see Lu, Darmofal (2004). 

Observe that the weight function (JJJ) is invariant with respect to rotations. Hence 
one might ask whether a result similar to (jHJ) holds for all symmetric weight functions. 
We conjecture that 

2d k 

(9) N min (2k + 1, d, 0) « — 

holds for all g and all k, hence the Moller bound is almost optimal. In this paper we 
prove this conjecture for k = 2 and k = 3, see Theorem 1 for more details. We also 
prove that the numbers N min (2k + 1, d, g) only mildly depend on the weight function 
g, see Theorem 2 for the details. 
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3 Some facts about the Smolyak method 



We study a special case of the Smolyak method, as we need it in the following. We 
also present methods with the upper bound (JBJ), since they are used (twice) for our 
new algorithm with the improved bound. We believe that this proof technique can 
be used to establish the conjecture © in full generality. Everything in this section 
is known or a minor modification of known results, see Novak, Ritter (1999). 

We construct cubature formulas to compute the integral (jSj) as follows. First we 
select quadrature formulas Uj, Uj, . . . to compute the one- dimensional integrals 

/ f(x) Qj(x)dx. 

These formulas should have the following properties: The formula C/j is exact for all 
univariate polynomials of degree mi, where 

(10) mi > 2i - 1. 

The formula £/j uses the knots Xj, the number rij = |Xj| of knots satisfies 

(11) ni<2i-l. 

We also assume that the are symmetric and "embedded" or "nested", i.e., 

(12) Xj" 1 C Xj for every i and j. 

By (lllip and (jllj) the weights of C/j are uniquely determined by its knots. Formulas 
with this property are often called interpolatory quadrature formulas. For simplicity 
we assume in this paper that the numbers rrii and rii do not depend on the coordinate 
j. The formula Uj, however, may depend on j. 

A product formula Ul 1 (g) • • • <S> U z d d needs . . . n« d function values, sampled on a 
grid. The Smolyak formulas A(q, d) are linear combinations of product formulas with 
the following key properties. Only products with a relatively small number of knots 
are used and the linear combination is chosen in such a way that the interpolation 
property for d = 1 is preserved for d > 1. The formula A(q, d) is defined by 

(13) A(q,d)= £ (-lr 111 -^,!,)-^ 1 ®---®^/), 

where q > d, i £ N d , and |i| = %\ + ■ — h 

The cubature formula A(q, d) is based on the sparse grid 

H(q,d)= U^x-xi;, 

|i|=9 

we use 

n = n(q, d) 
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to denote the cardinality of H(q, d). 1 In particular we have n(q, 1) = n q and we put 
n(0, 1) = no = 0. The recursion formula 

q-d+l 

(14) n(q + 1, d + 1) = n(q + 1 - s, d) ■ (n a - n s _ x ) 

8=1 

for n(q,d) is known, see Novak, Ritter (1999). 

Remark 1. Cubature formulas with high polynomial exactness are not often used 
if d is large, say d > 5. One major exception is the class of fully symmetric rules for 
the fully symmetric case, where also 

Qi = • • • = Qd- 

Fully symmetric cubature formulas were developed by Lyness (1965a, 1965b), Mc- 
Namee and Stenger (1967), Genz (1986), Cools and Haegemans (1994), Capstick 
and Keister (1996), Genz and Keister (1996) and other authors. The best results 
with respect to polynomial exactness are obtained by Genz (1986) and Genz and 
Keister (1996). The fully symmetric formulas from Genz (1986) and Genz and Keis- 
ter (1996) are of the Smolyak form ()13|) . Numerical integration with the Smolyak 
construction was already studied in Smolyak (1963). There are many other papers 
on the Smolyak method. The papers Gerstner, Griebel (1998), Novak, Ritter (1999), 
and Petras (2003) study the polynomial exactness of A(q,d). See also Novak, Rit- 
ter, Schmitt, Steinbauer (1999) and the recent survey on sparse grids by Bungartz, 
Griebel (2004). 

The following result is well known, see Corollary 1 of Novak, Ritter (1999). 

Lemma 1. Assume (JTUJ). Then A(d + k,d) has (at least) a degree £ = 2k + 1 of 
exactness. 

Now we present formulas for the number n(q, d) of knots that are used by A(q, d). 
We consider two cases, important for the following. 

The case rij = 2i — 1. 

Using ()14j) one obtains the recursion 

(15) n(q + l,d+l) = n(q, d + 1) + n(q, d) + n(q - 1, d) 

for q > d and n(q, 1) = 2q — 1 and n(d, d) — 1. Table 1 consists of numbers n(q, d) 
with minimal q such that n(q, d) > £, these numbers are called N(i, d). 

Using (fT3j) one can get an explicit formula for n(k+d, d), see Novak, Ritter (1999). 

1 Observe that some elements of the sparse grid might get a zero weight in the formula A(q, d). 
This would decrease the number of needed function values. Hence the "actual" number of needed 
function values for A(q,d) might be smaller than n(q,d). 
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Table 1: Number of knots for Smolyak's method with m = 2% — 1 



t 


AT ( K\ 

iv [Z, 5 J 


ATI 1 n"\ 

iv [Z, m ) 


AT l 1 K\ 

iv [Z, 10 J 


ATI o or\\ 
iv [Z, zll J 


AT 1 Ot^\ 
iV [I, ZD ) 


6 


i i 
1 1 


Zl 


ol 


41 




5 


61 


221 


481 


841 


1301 


7 


231 


1561 


4 991 


11521 


22151 


9 


681 


8 361 


39 041 


118 721 


283 401 


11 


1683 


36 365 


246 047 


982 729 


2 908 411 


13 


3 653 


134 245 


1 303 777 


6 814 249 


24 957661 


15 


7183 


433 905 


5 984 767 


40 754 369 


184327311 


17 


13 073 


1256 465 


24 331777 


214 828 609 


1 196 924 561 



Lemma 2. For every fceN and d G N we have 

Remark 2. Lemma 121 immediately implies 

(16) n(ife + d, d) < r + d ■ min(2 fc , 2 d ). 

The case = 2i — 1 for z 7^ 3 and n 3 = 3. 

If we take the Gaussian formulas Uj with 3 knots for Qj, then we already have 
exactness 5 and so we can take Uj = Uj and still have (jlOj) . Altogether we have 

(17) Hi = 2i — 1 for i ^ 3, n 3 = 3. 

Observe that in this case the sets X? are determined by the weights Qj, we cannot 
choose these sets. All the other sets Aj can be chosen arbitrarily for i > 2, but we 
still assume (fl2"j). Similarly as (fi3j) we now obtain from (fT4]) the recursion 

n(q + 2,d+l) = n(q + 1, d + 1) + n(g + 1, d) + n(q, d) 

-2n(q - 1, d) + An{q -2,d)- 2n(q - 3, d). 

With this simple modification we obtain the values of Table 2. 

Remark 3. Later the following will be important for the two versions of Smolyak's 
algorithm: In the case rij = 2i — 1 we can take arbitrary symmetric sets Aj, in 
particular we can take 

xl = ■ ■ ■ = xl 
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Table 2: Number of knots for method (|17|) 



f: 


AT 1 K\ 


at ( i n\ 
iv (£, 1U ) 


AT 1 1 K\ 
iv (£, to J 


AT l D on\ 
iv (£, zU J 


AT 1 OK\ 
iV [l, ZD ) 


o 
o 


it 


21 


ol 


41 


51 


5 


51 


201 


451 


801 


1251 


7 


151 


1201 


4151 


10 001 


19 751 


9 


401 


5 301 


27 701 


90 601 


227001 


11 


1003 


19 505 


146 507 


643 009 


2 040 011 


13 


2133 


63 805 


655 017 


3 775 769 


15 056 061 


15 


4 223 


188 745 


2 584167 


19111089 


94 680111 


17 


8113 


511625 


9 224 937 


85 920 449 


522 028 561 



We also can normalize the weights Qj in such a way that the U 2 have the form 

U*(f)=^f(-x)+p j f(0)+ 1 f(x), 

where 7 (and x) do not depend on j. In addition, we can choose the X 1 - in such a 
way that ||x|| 2 < a for each x e H(q, d), where a is the (given) radius of the domain 
Q of integration. This means that each rotation maps x to a point in Q. 

In the second case, however, we have to use the 3 GauB-knots for Xj = X?. 

Remark 4. Later we project the points H(q, d) of A(q, d) to a sphere of fixed 
radius. The origin is not projected. This projection reduces the number of points, 
the number of projected points n*(d + k, d) also depends on the sets X 1 -. We only 
need the second case, where n« = 2i — 1 for i 7^ 3 and n 3 = 3. In the case k = 2 and 
k = 3 one obtains 

n(d + 2,d) = 2d 2 + 1 and n*(d + 2, d) = 2d 2 

and 

n(d + 3,d) = (4d 3 -6d 2 + 20d + 3)/3 and n*(d + 3, d) = (Ad 3 - 6d 2 + 8d)/3. 

For the last formula observe that H(d + 3, d) contains 7 points of the form x = 
(a, 0, ... , 0) that are projected onto two different points, hence 

n(d + 3, d) = n*(d + 3,d) + 4d+ 1. 

It seems to be difficult to compute the smallest possible number n* (d + k, d) for 
general k, but it is clear that 

n(d + k,d) > n*(d + k,d) > 2 k 

Hence, for large d, we have n(d + k,d) ~ n*(d + k,d). 
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4 Known results for the Lebesgue measure 



Here we explain the best known upper bounds for N m[n (£, d, 1) that we found in the 
literature. Again we only discuss results for large d? 

The results for I e {3, 5, 7} are classical results that can be found in Stroud 
(1971): 

n = 2d for the degree £ = 3; this bound is sharp, A^" min (3, d, g) = 2d; 

n = 2d 2 + 1 for the degree £ = 5; 

n = (Ad 3 - 6d 2 + Ud + 3) /3 for the degree £ = 7. 

These results can be obtained with Smolyak's method, we explain the case £ = 7: 
First we take, as in ()17|). the values n 2 = n 3 = 3 and n 4 = 7. Now observe that 
the 4 new points of Xj are symmetric but otherwise arbitrary. Hence we can take 
(together with 0) the 5-point Gaufi rule with degree 9. This means that 2c? weights 
disappear and hence n is decreased by 2c? compared to the general situation of (JT7J). 

The best results (so far) for £ > 7 can be described in the following way: We 
use again the sequence m ; > 2i — 1 and so called "delayed Kronrod-Patterson- 
formulas". The rii are defined as follows: n\ — 1, n-x — n% — 3, = = = 7, 
n-j — • ■ ■ — n\2 = 15, 7ii3 = • ■ ■ = tt-24 = 31 and so on. Some of these numbers are 
larger than 2i — 1 and hence we can modify those n^, used by Petras (2003), to 

hi := min(nj, 2i — 1). 

In this way one obtains the values from Table 3, see Genz (1986) who obtained the 
same results. 



Table 3: Known values for the Lebesgue measure 



£ 


N(£,5) 


N{£, 10) 


N{£, 15) 


N(£, 20) 


N(£,25) 


3 


10 


20 


30 


40 


50 


5 


51 


201 


451 


801 


1251 


7 


141 


1 181 


4121 


9 961 


19 701 


9 


391 


5 281 


27671 


90 561 


226 951 


11 


903 


19105 


145 607 


641 409 


2037511 


13 


1733 


60 205 


642 417 


3 745 369 


14 996 061 


15 


3 263 


168 825 


2 473 287 


18 743 249 


93 755 311 


17 


5 983 


431265 


8 522 247 


82 703 329 


511676 911 



2 We illustrate this by an example. In the case d — 10 and I = 13 we will mention a method of 
Genz (1986) using n = 60 205 function values. In the same paper Genz presents another method 
using only n = 37 389. This method, however, uses more than 2 d points for general d and hence is 
not good for "large" d. 
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Remark 5. Observe that, up to now, there is nothing better known than the fully 
symmetric formulas that were introduced more than 40 years ago. We do not claim 
that the results of Table 3 are optimal for fully symmetric (or Smolyak) rules. It was 
proved by Petras (2003), however, that only minor improvements are possible if one 
uses Smolyak formulas. The same also holds for the more general fully symmetric 
formulas. For fixed i = 2k + 1 and large d, the number of points is (at least) of the 
order 

(18) jv(2fc + l,d,0)«l^, 

while the lower bound of Moller is only of the order 2 % r . Observe that (fTSJ) holds for 
all the versions of Smolyak's method that we presented here. 

Remark 6. By Lemma |2] we have the bound 

N(2k + 1, d, l)<[ k+ d d )- min(2 fc , 2 d ) 

for the Smolyak methods described here. For fixed k also Kuperberg (2004) obtains 
a bound of the form 

N(2k+i,d,i)< ( k+ /yc k . 

The constant C k is of the order 2 • k k ■ k\, much bigger than 2 k . However, Kuperberg 
(2004) obtains cubature formulas with positive (even equal) weights. This is a great 
advantage, in particular if the function values /(xj) are given only approximately. 
For a cubature formula Q n we define its condition number 



I ^ \ ||Qn||oo 2~2i= 



A cubature formula with positive weights has condition number cr(Q n ) = 1 if it is 
exact for the constant functions. The known Smolyak formulas of degree 5 and 7 
have a condition number of roughly d 2 and d 3 , respectively. See Remark 8 which 
also shows that our new formulas have roughly the same condition numbers. 



5 Cubature formulas for the sphere and for M^k 

In the following we need some known results for cubature formulas for the sphere. We 
use these results and the Smolyak method to construct efficient cubature formulas 
for the linear functional 

M d , k (f)= /(*) 

xeF(d,fc) 

where 

xGF(rf,fc) Xi G {±1, 0}, ^2x 2 = k. 



9 



Of course M^k itself is a cubature formula using 2 k (, ) function values, where k < d. 
The point is to find a cubature formula for M^k that is exact for polynomials from 
¥(2k + 1,0?) and uses only about 2(f) m 2d k /k\ points, which is the order of the 
lower bound of Moller. 

To achieve this we use two cubature formulas for the sphere that are exact for 
polynomials in F(2k + 1, d). The first formula is obtained from the Smolyak method 
for the Gaussian weight function (J7J) by projection onto the sphere of radius yk. It 
has the form 

(19) wM d>k (f) + Q r (f) 

where Q r (f) is a cubature formula with r = 0(G? fc_1 ) points and w > 0. In particular, 
we can take r < n* (d + k, d) — 2 k (^) with n* (d + k, d) from Remark |3J This leads to 

r < 2d for k = 2 and r < 2d 2 for k = 3. 

This works for any degree 2k + l of exactness. The second formula Q n {f) for k = 2, 3 
is taken from Mysovskikh (1968), see also Mysovskikh (1981). It uses 

n = d 2 + 3d + 2 points if k = 2 and d > 4 

and 

n = (d 3 + 9d 2 + Ud + 6) /3 points if k = 3 and d > 6. 

It follows that the formula w _1 (Q„(/) —Q r (f)) is a cubature formula for M d)k exact 
for polynomials from F(2k + 1, d) which uses at most 

(20) d 2 + 5d + 2 and (d 3 + 15d 2 + Ud + 6)/3 

points for k = 2 and k = 3, respectively. 

Let us finally explain how a Smolyak formula for the Gaussian weight function 
leads via projection onto the sphere of radius R = \fk to a cubature formula 

of the same degree of exactness. To this end, for r > 0, let u r be the surface measure 
on the sphere of radius r. Let also P be the radial projection from M d \ {0} onto 
jj>gd-i gj ven by p x _ / 1 1 x 1 1| . Furthermore, let 



QnCf) = X>/( X *) 



i=l 



be an arbitrary cubature formula which is centrally symmetric. Obviously, any 
Smolyak formula considered above has this property. We assume that Q n has degree 
of exactness 2k + 1 for the Gaussian weight function. Let x a = x" 1 . . . x a d d be a 
monomial of degree |a| = a.\ + . . . + a d = 2k. Using polar coordinates, we obtain 



x Q exp(— ||x|| 2 )(ix = / / x a <iu; r (x)e r dr 

'0 JtS^ 1 



POO P 

/ (r/R) d - 1+2k e- r2 dr / x Q ^(x) 

JO ii?S d -! 



c{R } d,k) \ x a c?u; fi (x). 
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We also have 

n n 
i=l i=l 

Whenever one of the points Xj = 0, we simply drop the corresponding term. Since 

Q„(x a )= / x Q exp(-||x||2)dx, 

we obtain that 

pg n (x a ) = / x Q rf^(x) 

where 

n 

PQ„(/) = $>/(P Xl ) 
i=l 

with 



1 R 2k c(R,d,k)' 

So PQn(f) is a cubature formula for the sphere RS^ 1 which is exact for homoge- 
neous polynomials of degree 2k. Since it inherits the central symmetry from Q n , it 
is also exact for homogeneous polynomials of degree 2k +1. Since any polynomial in 
F(2k + 1, d) restricted to RS d ~ l is a sum of two homogeneous polynomials of degree 
2k and 2k + 1, respectively, PQ n is exact for all such polynomials. 

If we choose the sets Xj in the construction of the Smolyak formula for the 
Gaussian measure equal, say Xj = X 2 = { — a, 0,a}, then the points x e aF(d, k) 
are present in the Smolyak formula and get equal positive weights. So the projection 
of this formula to the sphere RE> d ~ l has indeed the form (|19|). 

Remark 7. It will be important later on that the cubature formula derived for M^k 
uses only points on the same sphere of radius R = \fk where the points in F(d, k) 
live. 



6 Cubature formulas for general weight functions 

We now derive our main result which is formulated in the following theorem. 

Theorem 1. Let Q and g be as always and let k = 2,3. In the case k = 2 we assume 
d> 4, in the case k = 3 we assume d > 6. Then there exists a cubature formula Q n 
for l q d with degree 2k + 1 of exactness which uses at most 

(21) d 2 + 9d + l and (d 3 + 33d 2 + Ud + 3)/3 

points for k = 2 and k = 3, respectively. If the one- dimensional weight functions Qi 
are equal (the fully symmetric case) then the number of points can be reduced to 

(22) d 2 + 7d+l and (d 3 + 21d 2 + 20d + 3)/3 
for k = 2 and k = 3, respectively. 
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Proof. We start by describing how one can pass from the special cubature formulas 
for Md,k constructed in the preceding section to cubature formulas for general weight 
functions g as in the introduction. By proper scaling, we may assume that the radius 
of the domain O of integration is at least \fk. First, choose a Smolyak formula Q m 
for g that is exact for polynomials from F(2k + 1, d) and satisfies Xf = . . . = Xj = 
{ — 1, 0, 1}. Then Q m has the form 

(23) Q m = vM dyk + Q s 

for some v > and 

s = n(k + d,d) - 2 k 

In general, we have to use the case where rij = 2i — 1 for all % > 1. Then we obtain 

s = Ad + 1 and s = (18c/ 2 + 3)/3 

for k = 2 and k = 3, respectively Now we replace the part in (I23|) with the 
formula derived in the preceding section which uses at most as much points as given 
in ()20|) . By Remark [7| all points of the final cubature formula 

v ~ 

(Qn Qr) Qs 

W 

are in the interior of fl This cubature formula needs at most n + r + s function 
values. This leads to cubature formulas with 

d 2 + 9d + 3 and (d 3 + 33d 2 + Ud + 9)/3 

points for k = 2 and k = 3, respectively, which exceeds (|2Tj) by just two knots. 

A further reduction is possible if knots of Q n , Q r and/or Q s coincide. We explain 
how this leads to the reduced number of knots in (|22j) in the fully symmetric case. 
The reduction by two knots in the general case is achieved similar (and easier). 

To simplify notation, we denote by M r dk for r > the cubature formula 

MUf)= Yl /( x ) 

where 

x G F r (d, k) x, t E {±r, 0}, x 2 = kr 2 . 

Observe that M] fc = M dik . 

We further need some notation for fomulas derived from the simplex. Let 5* be 
a regular simplex with vertices in the unit sphere S d_1 . Let S dk be the cubature 
formula 

S" d \ k U)= E (/(x) + /(-x)) 

xeG r (d,k) 
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where G r (d, k) is the set of all projections of the centers of the (k — l)-dimensional 
faces of S onto the sphere of radius r. For the formulas of degree 7 we need one more 
cubature formula. Denote by Pij the (d + l)d points of the form 



1 3 

Pij = l v i + 1 v v 



where Vi and Vj are different vertices of the simplex. Then let H r (d) be the set of 
all rPy/IIPyll an d define the cubature formula by 

S r d (f)= £ (/(x) + /(-x)). 
xeH r (d) 

Finally, let Ud be the surface of S^ 1 . 

So assume now that Qi — ■ ■ ■ — Qd- We further assume without loss of generality 
that Q D [-1, l} d . We treat the de gree five and seven cases separately. 

Degree five. The projected Smolyak formula with degree of exactness 5 for the 
sphere S d_1 with d > 3 needs 2c? 2 points and has the form 

(24) ui Ml, + u 2 M X Jf 
with 

Ul = TdJdTTf* and U2 = WPif d - 

This formula can be found in Stroud (1971) or as formula 11) for the sphere in 
Mysovskikh (1981). 

The second formula with degree of exactness 5 for the sphere E> d ~ l with d > 4 
needs (d + l)(d + 2) points and has the form 

(25) V! Sl A + v 2 Sl 2 
with 

d(7-d) , 2(d-l) 2 

Vi = 2( d+ md+2f d and V2 = d{d+m d +2f d - 

This formula can be found in Mysovskikh (1968) or as formula 7) for the sphere in 
Mysovskikh (1981). 

Putting ()24j) and (|25|) together gives the following formula with degree of exact- 



ness 5 for M X J 2 2 : 

(26) -Ki+^4- U iMy. 

We also need a Smolyak type formula for the weight function g with degree of 
exactness 5 which has the form 

(27) a x M l Jf + a 2 M X Jf + a 3 Mj 1 + a 4 Q , 



13 



where Q (f) = f(0) and 7 e (0,1) \ {l/y/2}. The coefficients Oi, . . . , 04 can be 
derived either from the Smolyak construction or from direct computation using 
Sobolev's theorem which tells us that our formula has the required degree of exact- 
ness if it integrates the polynomials 1, X ^ - X ^ - X ^ X 2 correctly. This leads to a linear 
system of 4 equations for ai, . . . , 04 which has a unique solution. To minimize the 
number of knots we choose 7 = 1. 

1 /v^2 Anm ____ 

Finally, we replace M d 2 in formula (J27)) with the expression (J26|) . This leads to 
a formula 

(28) ai + 5^) + « 2 + a 3 M]/ 2 + « 4 Q 

which is exact of degree 5 for integration with respect to g with d > 4. The co- 
efficients cti, . . . , «4 can be directly derived using the polynomials 1, 
Alternatively, they are related to ai, . . . , 04 via 

2(d-l) 2 4-d 
(a + 1) 2 2 

Observe that we have chosen our formulas so that the final number of knots is 
d 2 + Id + 3. This can be further reduced to 

d 2 + 7d + 1 

if we choose one of the vertices of the regular simplex S as the unit vector (1, 0, . . . , 0). 
Observe also that in the case d = 7 the number of knots reduces even further. 

Degree seven. Let us now derive a formula with degree of exactness 7, i.e., k = 3. 
The projected Smolyak formula with degree of exactness 7 for the sphere S d_1 with 
d > 3 needs (4<i 3 — 6d 2 + 8d)/3 points and has the form 

(29) « a Ml + u 2 M l Jf + u 3 M]!f. 

This formula can be found in Stroud (1971) or as formula 21) for the sphere in 
Mysovskikh (1981). 

The second formula with degree of exactness 7 for the sphere E> d ~ l with d > 6 
needs (d 3 + 9d 2 + 14d + 6)/3 points and has the form 

(30) v t Sj A + v 2 S\ 2 + v 3 Sj 3 + v 4 Sl 

This formula can be found in Mysovskikh (1968) or as formula 13) for the sphere in 
Mysovskikh (1981). 

Putting (|2U|) and (j3~Uj) together gives the following formula with degree of exact- 
ness 7 for M^^: 



(31) - U S l d>1 + v 2 Sl 2 + v 3 S l dt3 + v 4 S l d - Ul Ml, - u 2 M l Jf 

u 3 \ 
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We also need a Smolyak type formula for the weight function q with degree of 
exactness 7 which has the form 

(32) ai M]'f + a 2 M l Jf + a 3 M x Jf + a 4 M% + a 5 Mj\ + a 6 Af£ + a 7 Q , 

where Qo{f) = /(0) and the numbers 71 and 72 and l/v3 are pairwise different, 
between and 1. To minimize the number of knots in the following we choose 
71 = l/y/2 and 72 = 1. 

1 / \Zs _____ _____ 

Finally, we replace M d 3 in formula (|3*2*j) with the expression (|31j). This leads to 

a formula of the form 

«i («i Sji + V2 S} }2 + v 3 Sl >3 + v A Sjj + 

The constants «i, . . . , a 7 can be determined by using the 7 polynomials 1, xf, xf, 
2 and xj. Observe that we have chosen our formulas so that the 
number of knots is 

(d 3 + 21d 2 + 20d + 9)/3. 
This can be further reduced to 

(ci 3 + 2Li 2 + 20d + 3)/3 

if we choose one of the vertices of the regular simplex S as the unit vector (1,0,. ..,0). 

□ 

Table |U contains the number of function values for fully symmetric weight func- 
tions. Observe that for £ = 7 we have to assume d > 6. 

Table 4: New values for fully symmetric weight functions 

£ N(£,5) N(£,10) N(£,15) N(£,20) N(£,25) N(£,50) N(£,W0) 
5 61 171 331 541 801 2 851 10 701 

7 - 1 101 2 801 5 601 9 751 59 501 404 001 



It is interesting to compare these values with the lower bound (0) of Moller, see 
Table El 

Table 5: Moller's lower bound 

£ N(£,5) N(£,10) N(£,15) N(£,20) N(£,25) N(£,50) N(£, 100) 
5 31 111 241 421 651 2 551 10101 

7 80 460 1390 3120 5 900 44 300 343 600 
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Remark 8. For the cube [— 1, l] d with Lebesgue measure, Tables and [7| contain 
the coefficients a; and on in the cubature formulas (|2*7j) . (J2*8|l . The values of 

■Ui, . . . , V4 and Mi, «2, ^3 for the degree 7 formula can be found in Mysovskikh (1981). 



Table 6: Coefficients for the degree 5 formulas (|27|) and (|28j) 





2 




3 






4 


22 


2d 




1 




2d 2 


37c2 


45 


9 




30 




9 


45 


<i 


17 


22 




2d 


2d 2 


37d 


18 


90 


45 




9 


9 


45 



+ 1 
+ 1 



Table 7: Coefficients for the degree 7 formula 



1 _7_ 6? £0 _ 

8 20 4 70 20 ~ t ~ 4 45 



2 d di gg 



16d 
45 



21 



£(_ , 5cT _ t>59rf _i_ i 
' 6 ~ l ~ 9 630 



Remark 9. Victoir (2004) and Kuperberg (2004) describe, in particular, methods 
for I = 5 and positive weights. For d = 100 Victoir has n = 4 12 = 16 777216 and 
this was further improved by Kuperberg to n = 65 536 points with positive weights. 
See the discussion in Kuperberg (2004). 

For general weights the old record was 20 001, see (j3J). Our method needs 10 701 
function values, the lower bound of Moller is 10 101. 

7 Independence of the weight function 

We now use the Smolyak formulas to show that, for any fixed k, the minimal number 
of knots needed by a cubature formula of degree 2k + 1 does not essentially depend 
on the weight function. Since the Moller lower bound is of order d h , the following 
theorem shows that the difference can only be in the lower order terms. 

Theorem 2. Let 12 and qv\ j = 1,2, be two regions and weight functions in M d 
as described in the introduction. For k = 2, 3, . . define 
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Then 

\N min (2k + 1, d, g {1) ) - N min (2k + 1, d, g {2) )\ < c k d k ~ l 

for all d > k. 

Proof. Without loss of generality, we assume that the cube [— 1, l] d is contained 
in the interior of f^ 1 -* and Q( 2 \ We choose a cubature formula Q n for g^ exact 
for polynomials in F(2k + 1, d) with n = N m - m (2k + 1, d, g^ 1 '). By proper scaling if 
necessary we may now assume that the knots of Q n are in the interior of Q( 2 \ We 
also choose, for j = 1,2, Smolyak formulas 

Q s ™ l = w 3 M d , k + Q r] 

for gV' of degree 2k + 1 with wj > 0. To assure their existence, we have to work 
with the case rti = 2i — 1 for all i. In this case we can also arrange that the knots of 
Q r are contained in [—1, l] d . Then, for d> k, the estimate 

(33) r,<-*( d t k )- 2 t 
follows from JTSJ). Now 

-Q ri )+Q 

Wi 

defines a cubature rule for g^ exact for polynomials in P(2fc + 1, d) with at most 
n + r\ + r 2 knots. Observe that all the knots used are in the interior of fi (2) . By (1551) . 
to prove the theorem it is enough to verify the elementary inequality 

2*( d+k \-2*( d \<J^ i ^ 



k J \kj ~ (k- 1)! 

for d > k, which is equivalent to 

(d + k){d + k - 1) . . . (d + 1) - d(d - 1) . . . (d - k + 1) < k 2 k d k ~\ 

Since the left-hand side of this inequality does not exceed (d + k) k — (d — k) k , this 
is an immediate consequence of 



( d + - (d - A;) fc = 2 ^ ^ d fc " J F < 2d fc - 1 A; ^ ^ = fc 2 k d 



fc jfe-i 



0<i</c x ' 0<i<fe 

odd i odd 



□ 

Remark 10. Similarly, it can be shown that 

|iV min (2£; + 1, d, (i d ) - A^ min (2A; + l,d,g)\< c k d k -\ 

where fid is the surface measure on the sphere § d_1 and g is a weight function as in 
Theorem EJ 

Acknowledgment. We thank two anonymous referees for helpful comments. 
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